Practality of the non-equilibrium stationary states of open 
volume-preserving systems: I. Tagged particle diffusion 
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Deterministic diffusive systems such as the periodic Lorentz gas, multi-baker map, as well as spa- 
tially periodic systems of interacting particles, have non-equilibrium stationary states with fractal 
properties when put in contact with particle reservoirs at their boundaries. We study the macro- 
scopic limits of these systems and establish a correspondence between the thermodynamics of the 
macroscopic diffusion process and the fractality of the stationary states that characterize the phase- 
space statistics. In particular the entropy production rate is recovered from first principles using 
a formalism due to Gaspard [J. Stat. Phys. 88, 1215 (1997)]. This article is the first of two; the 
second article considers the infiuence of a uniform external field on such systems. 

PACS numbers: 05.45.-a,05.70.Ln,05.60.-k 



I. INTRODUCTION 

The assumption that the microscopic dynamics of me- 
chanical systems obeying Newton's equations are mixing 
offers a mechanism by which the entropy can increase to- 
ward its equilibrium value, as described by Gibbs in 1902 
[l[ . The mixing would indeed allow coarse-grained proba- 
bilities to reach their equilibrium values after a long time, 
a result that has received a rigorous meaning in the con- 
text of modern ergodic theory 0] . The use of the coarse- 
grained entropy of a physical system, which is the second 
ingredient of Gibbs' mechanism, is justified by the fact 
that, if the entropy should be given according to Boltz- 
mann by the logarithm of the number of complexions of 
a system, then the introduction of cells of non-vanishing 
sizes is required to perform the counting of complexions 
in systems described by continuous coordinates. 

The program set up by Gibbs has taken on a new 
perspective in recent years with its systematic applica- 
tion to chaotic, deterministic volume-preserving dynam- 
ical systems sustaining a transport process of diffusion, 
which therefore fall under the scope of Liouville's the- 
orem [3, i, i, i, 0, i, II- In previous works, both 
non-equilibrium stationary states and relaxation to equi- 
librium were considered, with the common thread that 
transport processes are intimately connected to singular- 
ities of the non-equilibrium measures, characterized by 
phase-space distributions with fractal properties. 

A specific example of the systems of interest is the 
periodic Lorentz gas, in which moving particles diffuse 
through a lattice interacting only with fixed scatterers. 
A useful simplification of this process is the multi-baker 
map, which is a simple model of deterministic diffusion 
for a tracer particle with chaotic dynamics. Other higher 
dimensional examples are spatially periodic extensions 
of systems with many interacting particles, where the 
motion of a tagged particle is followed as it undergoes 



diffusion among the cells. 

In refs. H, Q , these models were studied in the con- 
text of relaxation to equilibrium. It was shown in these 
papers that the initial non-equilibrium distribution func- 
tion rapidly develops a fractal structure in phase space 
due to the chaotic nature of the dynamics. This stuc- 
ture is such that variations of the distribution function 
on arbitrarily fine scales grow as the system evolves in 
time. The final stages of the approach to equilibrium are 
then controlled by the decay of fractal, microscopic hy- 
drodynamic modes of the system -in this case diffusive 
modes- which decay with time as exp(— 2?fc^t), where k 
is a wave number associated to the macroscopic hydro- 
dynamic decay, T) is the diffusion coefficient, and t is the 
time. It is possible to express the rate of entropy pro- 
duction in this final stage in terms of measures which are 
determined by the non-equilibrium phase-space distribu- 
tion, specified by the fractal hydrodynamic modes. The 
main result is that one obtains by this method exactly the 
expression for the rate of entropy production as given by 
irreversible thermodynamics for these systems [lOj . The 
source of this agreement can be traced to the role played 
by the fractal hydrodynamic modes, both for requiring 
a coarse graining of the phase space to properly incor- 
porate the effects of their fractal properties on entropy 
production in the system, as well as for describing the 
slowest decay of the system as it relaxes to equilibrium. 

The purpose of this paper is to consider in further de- 
tails the non-equilibrium stationary states of multi-baker 
maps and Lorentz gases which occur when their bound- 
aries are put in contact with particle reservoirs, as well 
as to extend these considerations to tagged particle diffu- 
sion in spatially periodic systems of interacting particles. 
In the macroscopic limit, we establish the connections 
between the statistics of these deterministic systems and 
the phenomenological prescriptions of thermodynamic. 
Though the stationary states of multi-baker maps and 
Lorentz gases have been considered in some details, see 
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especially [Tl|, the problem of computing the entropy 
production associated to the non-equilibrium stationary 
state has thus far been limited to the example of the 
multi-baker map Here we emphasize the similari- 
ties between the stationary states of multi-baker maps, 
Lorentz gases, and spatially periodic many particle sys- 
tems, undergoing steady mass flows and show how the 
formalism described in [9i] yields an ab initio derivation 
of the entropy production rate in all these systems. 

This paper is the first of two. In the second paper 
[l^ . we will consider the Galton board [H, which 
is to be understood in this context as a periodic finite 
horizon Lorentz gas with a uniform external field and 
no dissipation mechanism. As was recently proved by 
Chernov and Dolgopyat [3, [3] , this system is recurrent 
and therefore has no drift. We will show how one can 
characterize equilibrium and non-equilibrium stationary 
states of such a system, much in the same way as with 
the periodic Lorentz gas we consider in this paper. 

The plan of the paper is as follows. In Sec. |TT1 we re- 
view the phenomenology of non-equilibrium diffusive sys- 
tems and their entropy production. Then, starting from 
the Liouvillian evolution of phase-space distributions, we 
construct the non-equilibrium stationary states of open 
periodic Lorentz gases in Sec. lIIIl as well as of multi-baker 
maps in Sec. IIVI and show how such stationary states 
naturally separate between regular and singular parts. 
While the regular parts have the form of a local equilib- 
rium and allow, in the continuum limit, to retrieve the 
phenomenological solution of a diffusive system undergo- 
ing a mass flow, the singular parts have no macroscopic 
counterpart and encode the dynamical details of the sys- 
tems' diffusive properties. In Sec. |Vl these results are 
extended to systems of many interacting particles with 
non-equilibrium flux boundary conditions. Under the as- 
sumption that the dynamics of these systems is chaotic, 
their stationary states are shown to have properties sim- 
ilar to those of multi-baker maps and Lorentz gases. The 
singular part of the stationary states is at the origin of 
the entropy production rate, as is explained in Sec. IVIl 
where the phenomenological entropy production rate is 
retrieved from ab initio considerations. Conclusions are 
drawn in Sec. IVIII 



II. PHENOMENOLOGY 

The thermodynamics of deterministic models of dif- 
fusion such as the Lorentz gas and multi-baker map to 
be discussed below was extensively reviewed by Gaspard 
et. al in [17]. In these models, independent tracer par- 
ticles mimic matter exchange processes of binary mix- 
tures, i.e. the process of mutual diffusion between the 
light tracer particles and the infinitely heavy particles 
that constitute the background. 

In such systems, the irreversible production of entropy 
arises from gradients in the density of tracer particles. 
For a dilute system, the resulting thermodynamic force 
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where V{X,t) denotes the density of tracer particles at 
position X and time t. Here the Boltzmann constant is 
set to unity. Close to equilibrium, in the linear range of 
irreversible processes, the corresponding current is pro- 
portional to the thermodynamic force and given by Pick's 
law of diffusion according to 



J{X,t) = -VVV{X,t), 



(2) 



where V is the coefficient of diffusion of tracer particles 
(assumed to be uniform). In this limit, the condition of 
mass conservation, dtP{X, t)+^-J{X, t) = 0, transposes 
into the the Fokker-Planck equation. 



dtV{X,t) ^VV^V{X,t). 
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The product of the thermodynamic force and current 
yields the rate of local irreversible production of entropy. 
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Our goal in the next sections will be to analyze the 
statistics of deterministic models of diffusion, identify the 
conditions under which their statistical evolution reduces 
to the evolution of macroscopic densities as described by 
Eq. ([3]) and analyze the microscopic origins of the pro- 
duction of entropy prescribed according to Eq. (j4|) , which 
manifest themselves in the fractality of the stationary 
states of these models. 



III. OPEN PERIODIC LORENTZ GAS 

In this section, we consider a two-dimensional finite 
horizon periodic Lorentz gas with hexagonal symmetry, 
in the shape of a cylindrical Lorentz channel of length L, 
with periodic boundary conditions in the y-direction and 
absorbing boundaries at a; = ±L/2. Though some of the 
results presented in this section are original, for the most 
part this material is a review of existing results that can 
be found in Ref. [ll[ and references therein. 

The geometry of the channel is displayed in Fig. [T] with 
a typical trajectory. The channel consists of a cylinder 
of length L = Nl and height ^/il with disks ]D)„, -N < 
n < iV, of radii a, '\/3/4 < a /I < 1/2. Assuming N even, 
the disks' centers take positions 



(n//2,0), nodd, 
n even, 



(5) 




FIG. 1: Cylindrical Lorentz channel with periodic boundary conditions at the upper and lower borders. The length of the channel 
is L, with A'^ identical cells of widths I and heights V31/2, Nl = L, each containing two disks. The disks have radii a — 0.44/, and 
are positioned at points {x,y) = (nZ/2,0), n = -A''+l, -A^ + 3, . . . ,iV- 1, and {nl/2, ±V3l/2), n = -N, -N + 2, ...,N-2,N. 
The trajectory displayed is that of a particle which initially starts at unit velocity from the left-hand boundary, and is eventually 
absorbed as it reaches the right-hand boundary. Arrows on the upper and lower borders indicate points where the trajectory 
winds around the cylinder. 



where the disks at y — ±^/3l/2 are identified. We will 
also denote by I„ the cylinder region around each disk 



to the outgoing normal to the disk after the collision. 



I„ = {{x,y) I {n-l/2)l/2<x<{n+ 1/2)1/2}. (6) 

Thus the interior of the cylinder, where particles propa- 
gate freely is made up of the union U^^_^I„ \ ID)„. 

The associated phase space, defined on a constant en- 
ergy shell, is C = U^^_^C„, where C„ = (g) [I„ \ B„] 
and the unit circle §^ represents all the possible velocity 
directions. Particles are reflected with elastic collision 
rules on the border dC, except at the external borders, 
corresponding to x = ±L/2, where they get absorbed. 
Points in phase space are denoted by F, F = (x, y, Vx, Vy), 
with fixed energy E = {v^ -\- Wy)/2, and trajectories by 
$*r, with <i>* the flow associated to the dynamics of the 
Lorentz channel. 

The collision map takes the point F — {x,y,Vx,Vy) € 
dC to $'^F — {x' ,y' ,v'^,Vy) e dC, where r is the time 
that separates the two successive collisions with the bor- 
der of the Lorentz channel 3C, and {Vx,Vy) is obtained 
from {vx,Vy) by the usual rules of specular collisions. 
Given that the energy E is flxed, the collision map op- 
erates on a two-dimensional surface, which, given that 
the collision takes place on disk n, is conveniently pa- 
rameterized by the Birkhoff coordinates (</>„, where 
(pn specifies the angle along the border of disk n that the 
trajectory makes at collision with it, and ^„ is the sinus 
of the angle that the particle velocity makes with respect 
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A. Non-Equilibrium Stationary State 

In order to set up a non-trivial stationary state, we 
assume that a flux of trajectories is continuously flow- 
ing through the boundaries. This can be achieved by 
putting the boundaries in contact with particle reservoirs 
such that the phase-space density at x = zLL/2 is flxed 
according to constant values. 



p(r,t) 



= ±L/2 



P±, 



(7) 



corresponding to particle injection rates from the two 
ends of the channel at cc = ztL/2. All the injected par- 
ticles have the same energy E — (u^ + w^)/2 and uni- 
form distributions of velocity angles. The evolution of 
the phase-space density p is otherwise specified by Li- 
ouville's equation, or, equivalently, by the action of the 
Frobenius-Perron operator, P*, determined according to 



p(F,t) = PV(r,0)= / dF'<5(F-<i>*F')p(F',0), 
Jc 



(8) 



A remarkable result [l8| is that the invariant solution of 
this equation, compatible with the boundary conditions 
([7]), is given, for almost every phase point F, by 



Pin = 



P++P'- , p+- p- 
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In these expressions, a;(r) denotes the projection of the phase point F on the horizontal axis, and T(r) is the time 
it takes the phase point F to reach the system boundary backward in time, i.e. such that a;(<i>~^*^'"^F) — ±L/2. As 
written in the second hne, the time T(F) elapses in K{T) successive collisions separated by time intervals r„, such 
that tk is the time elapsed after k collision events, viz. X]j=i ~ ^k- Obviously ^^(r) — T(T). The difference 
a;($~*'=F) — x{^~^''-^r) is the displacement along the x axis between two successive collisions. 

In the limit of infinite number of cells N, the number of collisions for the trajectory to reach the boundaries becomes 
infinite so that the invariant density can be written 

P(r) = + |:r(F) + f;[x(<i>-*^F) - :r(<i>-*-^F)]| . (10) 



The two contributions to this expression represent the 
mean linear density profile along the axis, plus fiuctu- 
ations. Notice that the linear profile is itself the sum 
of two contributions, the first one being the equilibrium 
density and the second one a gradient term. These sta- 
tionary states are known after Lcbowitz and MacLennan 
As noted in [l^, the fiuctuations form the sin- 
gular part of the invariant density, which builds up ac- 
cording to which of the two boundaries the phase point 
is mapped to. 



B. Continuum Limit 

The linear part of the invariant density profile Q has 
its origin in the diffusion process which takes place at the 
phenomenological level. Indeed, in the continuum limit, 
the phase-space density reduces to the macroscopic den- 
sity distribution V{X, t), whose evolution is described by 
Eq. ([3]), with macroscopic position variable X and diffu- 
sion coefficient, V, which can, in principle, be determined 
from the underlying dynamics. 

In order to obtain Eq. ([3]) from the Liouvillian evolu- 
tion, Eq. (O, we let ; ^ and iV -> oo with Nl = L 
fixed. All other quantities are fixed, in particular the 
particle velocity, which is taken to be unity in the length 
units of L. In that limit, the macroscopic density can be 
obtained from the phase-space density by averaging over 
phase-space regions such that the position variables x(F) 
are identified with the macroscopic position variable X, 
here denoted Xn and identified with the cell d, 

V{X„,t)^j [ dVp{T,t). (11) 

Note that the volume element is here and throughout 
assumed to be properly normalized, i.e. dT = I. 

Given non-equilibrium fiux boundary conditions with 
the stationary state (|10p . the only contribution to the 
particle density V{Xn) is the linear part, which, in the 
continuum limit, yields the stationary solution to Eq. ([3]), 
namely 

V{X) = '^y- + -(n - V-) , (12) 



where 7-'± is equal to p± up to a volume factor, and spec- 
ifies the boundary conditions 7^(±L/2). 

The singular part of the invariant density, on the other 
hand, has no macroscopic counterpart. Its integration 
over the cell indeed vanishes by isotropy of the underlying 
dynamics. Nonetheless, as we will demonstrate shortly, it 
bears essential information pertaining to the microscopic 
origin of entropy production. 

In that respect, the phenomenological entropy produc- 
tion rate ([?]) can be written 

dt " ^ v{x) ' ^^-^^ 

where we denoted T'cq — [V^ + V-)/2. 



C. Numerical Results 

We can verify numerically that the invariant density of 
the Lorentz channel has the form ([9]), with a linear part 
given by the stationary solution of the Fokker-Planck 
equation (fT^ . and a singular fluctuating part. For the 
sake of the numerical computation, we let L = 1 and 
consider a channel with 2N + 1 disks, N = 25, the left- 
and right-most ones being only half disks. The macro- 
scopic position X, —L/2 < X < L/2, will be identified 
with the cells C„ about the corresponding disks, letting 
Xn = nl/2, —TV < n < N. We then fix the injection rates 
p± and let p_ > and /o+ = so as to have V- = 1 and 
V+ ~ for the corresponding macroscopic quantities. 

Having fixed the geometry of the channel and the in- 
jection rates, we compute the invariant phase-space den- 
sity p(F) in terms of the corresponding density of the 
collision map, also called the Birkhoff map, which takes 
trajectories from one collision event with a disk to the 
next one. This is indeed much easier since the numerical 
integration of the Lorentz gas relies on an event driven 
algorithm corresponding to integrating the collision map. 
The conversion between the two phase-space densities is 
here trivial because the time scales are uniform. 
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Thus the Hnear part of the invariant density is com- 
puted by considering a large set of initial conditions in- 
jected into the channel from the left end and integrating 
them until they reach either ends of the channel, thus ex- 
iting the system. In the meantime, we record the num- 
ber of collisions each particle performs with each disk, 
the ensemble average of which approximates the station- 
ary distribution V{Xn). Provided N is sufficiently large 
and n is not too close to the channel boundaries, this 
ensemble average is expected to approximate V{X), as 
specified by Eq. p^ . The result is displayed in Fig. [2] 
and was adjusted by an overall constant so as to match 
V- = 1 (and V+ = 0). The agreement is very good. 




FIG. 2: Non-equilibrium stationary density of the Lorentz 
cliannel obtained for a cylinder of length L = 1, with 51 disks 
(iV = 25). The solid line is V{X) = 1/2 - X, the solution 
of Eq. (Ull), with = 1 and P+ = 0. Notice that the left- 
and right-most disks are only half disks. Hence the strong 
boundary effects. 

The singular part of the density can also be computed 
in terms of the statistics of the collision map. For the 
sake of characterizing this quantity, we notice that the 
equilibrium Lorentz gas (a single disk on an hexagonal 
cell with periodic boundary conditions) preserves the Li- 
ouville measure 

dT = dx dy dvx dvy , 

= dEdtd(j)dC, (14) 

where (j> and ^ are the Birkhoff coordinates. 

For the open Lorentz gas, we associate to each disk 
a pair of Birkhoff coordinates, (0n,Cn), and record each 
collision event in a histogram. The results, displayed in 
Fig. [3J show the fractality of the fluctuating part of the 
invariant density, with boundary effects vanishing expo- 
nentially fast with respect to the distance of the disk to 
the channel boundaries. 

We will see in Sec. |VT] that the fractality of the sta- 
tionary measure displayed in Fig. [3] is responsible for the 
entropy production rate that is associated to the mass 
current in the open Lorentz gas with flux boundary con- 
ditions. 



IV. MULTI-BAKER MAP 

The analytic derivation of the entropy production rate 
(|13p . relying on the fractality of the invariant measure 
pn)) . can be achieved with the multi-baker map. This 
is a simple model of deterministic diffusion that can be 
thought of as a caricature of the dynamics of the Lorentz 
gas described above. This model was originally intro- 
duced in [2l| , its non-equilibrium stationary state and re- 
lation to thermodynamics analyzed in \§\ , and the deriva- 
tion of the entropy production rate described in 

The basic idea which underlies the similarity with the 
Lorentz channel is threefold. First, let the positions take 
discrete lattice positions, like the cells C„ of the chan- 
nel; second, assume that the collision times t„ are all 
identical and denoted by t, which we leave arbitrary for 
now; and, third, replace the reflection rules at the col- 
lisions by a simple Bernoulli-type angle-doubling rule, 
which decides which direction the particle goes to. As- 
suming the system length L to be an integer multiple 
of the cell size I, which for the sake of the argument we 
take to be L = {N + 1)1 with N even, tracer dynamics 
of the Lorentz channel are replaced by the mapping on 
{-N/2, . . . , N/2} ® [0, l]^, given by 

Bo : {n,x,y) i-^ (15) 

f (n - l,2x,y/2) , 0<x <l/2, 

\ln + l,2x~l/{y + l)/2), l/2<x<l, 

with the convention that trajectories end when mapped 
to n outside of {—N/2, . . . , The map is displayed 

in Fig. m We point out that the multi-baker map has 
a time-reversal symmetry with respect to the operator 
S : (n, X, y) ^ {n,l — y,l — x), namely S o Bq — Bq^ o S. 
The map (fT5)) includes a displacement from cell to cell, 
which we denote by A, 

Hn,.,y)^[-\l lt<x<l: (16) 

Because Bq is a Bernoulli map, a point F = (n, x, y) 
in phase-space may be thought of as coding an infinite 
sequence of equi-probable nearest-neighbors jumps both 
in the past and future, with initial condition at position 
nl. Moreover the variables {x,y), play here the role of 
the Birkhoff coordinates ((/>, ^) of the Lorentz gas. 



A. Non-Equilibrium Stationary State 

Given flux boundary conditions, a non-equilibrium sta- 
tionary state similar to that of Eq. ^ sets in, with a 
singular part given in terms of the sum of the succes- 
sive displacements the tracer makes before it reaches the 
outer boundaries and gets absorbed. The x position now 
becomes a lattice coordinate n and the change in posi- 
tions between collisions the displacement ((T5)) . so that 
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FIG. 3: Fractal structure of the phase portraits of an open cylindrical Lorentz channel with absorbing boundaries at its ends. 
The plots are histograms computed over grids of 500 x 500 cells, and counting the average collision rates of many trajectories 
in every cell. Disk 50 is the one before last. The geometry of the channel is that shown Fig. [T] Here particles are injected at 
the left boundary only. The color white, associated to higher counting rates, corresponds to the injection of particles from the 
left boundary. Black areas, on the other hand, correspond to absorption at the right boundary. Hues of gray are associated 
to phase-space regions whose points are mapped backward to both left and right borders. The corresponding overall densities 
([TT]) are shown in Fig. [S] 
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n - 2 ri - I 



II + 1 11 + 2 



FIG. 4: (Color online) Multi-baker map ([T5|l . Each cell has 
area l^, with coordinates {x,y), and is labeled by an integer 
-N/2 <n< N/2. 



the stationary state density reads 

p(r) 

P+ + p- , - p- 



(17) 



Kir) 

(r)/ + ^ A(So-^T) 

k=l 



where n(T) denotes the projection of F along the integer 
axis. 

Since the sum over the successive displacements is a 
highly singular function of F, it is convenient to introduce 
a better behaved characterization of the invariant density 
in terms of cumulative functions. Thus let F = {n,x,y) 
with <x,y <l. We define the set(29t 

C„(x, y) = {(x', y')\0<x' <x,0<y' < y}. (18) 

For X — y — I, the cylinder set corresponds to the whole 
cell, C„ = C„(^, / ). The cumulative measure of this set is 
defined as 



Pn{x,y) = 



C„(x,y) 
P++P- 

2 

P+ - P- 
L 



dTpiV) , 



nl 

{P+-P-)J^ 

K(T) 



/ dF ^ A(So-^T) 

JC„{x,y) 



Here we again assumed normalization of the volume ele- 
ment, dr — dxdy/P. 

It can be shown [3|] that the expression above reduces 

to 



xy 



Pnix,y) ^ Pn-pT + 



{P+- P-) 
L 



xT„ (I) 



where the linear part of the invariant density is 



1 nl 
2{P+ + P-) + {P+ - P-)^ ^ 



(20) 



(21) 



and T„ is the incomplete Takagi function [3|], which can 
be defined through the functional equation, here with 

ye [0,1], 



Tn{y) 



2/+ir„+i(2y), 0<2/<l/2, 
l-y + ir„_i(2y-l), l/2<y<l. 



This equation can be solved recursively using dyadic ex- 
pansions of X and the boundary conditions T^i^/2{x) — 0. 
Notice that the size I of the jumps A was extracted from 
the Takagi functions and absorbed into x {= I x x/l) . 
We remark that verifies 



1 



1 



Pn — + -^Pn+l, 



(23) 



which describes the invariant statistics of a uniform ran- 
dom walk with jump probabilities 1/2 to nearest neigh- 
bors. Given the length / of the jumps, and their rate r, 
the diffusion coefficient for this process is D = P/(2t). 
The proper diffusive scaling is recovered provided t ^ P 
when r, / — !■ in the continuum limit. 

In the continuum limit, the number of steps that sep- 
arate phase points from the boundaries becomes infinite 
so that the invariant cumulative measure can be written 
in terms of the (complete) Takagi function. 



/ N xy [p+- p_) fy 

Pn[X,y) = Hn-j2 H -j^ xT' 



(24) 

The Takagi function [2^, displayed on the top panel of 
Fig. O is the solution of the functional equation 



Tiy) 



y + hT{2y), 



0<y < 1/2, 



l-y+^T{2y-l), 1/2 < y < 1. 



(25) 



(22) 



This equation, due to de Rham '23'|, is one of many 
equivalent representations of the Takagi function, T{y) = 
J2n=o 2"" |2"y - [2"y -I- 1/2] I, where [x] stands for the 
maximum integer not exceeding x. It is an everywhere 
continuous function, but has no bounded derivative. 

The singularity of the fluctuating part of the invari- 
ant density (jl7p with respect to the y component can 
be visualized by measuring the local slopes of the Takagi 
function. The dependence on the x component on the 
other hand is trivial. The resulting histogram, displayed 
on the bottom panel of Fig. [5l can be compared to those 
shown in Fig.[3]for the Lorentz channel. The main differ- 
ence is that while the fractal structures are curved in the 
open Lorentz system, they are parallel to the x-direction 
in the multi-baker map. The reason is that the fractal 
structures are smooth with respect to the unstable direc- 
tions, which, for the multi-baker map, are always directed 
along the x-axis, while their direction varies from point 
to point in the Lorentz system. 



B. Continuum Limit 

As with the Lorentz gas, we let / — > and N — > oo, in 
the continuum limit, keeping L = {N + 1)1 constant. In 
this limit, the macroscopic particle density is to be identi- 
fied with /i„, V{Xn — nl) = On the scale of /i„, the 
singular fluctuations embodied by the Takagi function 
disappear. However, as argued in these fluctuations 
are responsible for the positiveness of the entropy produc- 
tion rate, given according to phenomenology by Eq. (fie 
We will comeback to this identification in Sec. IVII 



8 




1 1 3 1 

4 2 4 



y 




Q L L 11 I 

4 2 4 

y 



FIG. 5: (Top) The Takagi function, Eq. ([25}, characterizes 
the fluctuations of the non-equilibrium stationary state as- 
sociated to the multi-baker map under flux boundary condi- 
tions. (Bottom) The singularity of the invariant density of 
the multi-baker map (|17|l can be visualized by measuring the 
local slopes of the Takagi function (|25|) at a given scale, here 
dy — 1/512. This histogram can be compared to those ob- 
tained for the Lorentz channel. Fig. [3] 



V. SPATIALLY PERIODIC SYSTEMS OF 
INTERACTING PARTICLES 

The results presented in Sees. IIIHIIVI for diffusive sys- 
tems of non-interacting tracer particles can be extended 
to spatially periodic systems of many interacting parti- 
cles where tagged particles -assumed to be independent 
of each other- undergo a process of self-diffusion. The ar- 
guments below provide an extension to non-equilibrium 
stationary states of the results presented in [9;] for tagged 
particle diffusion in interacting particle systems under- 
going a relaxation to equilibrium. 

We follow the diffusion of tracer particles in a periodic 



array of cells, each containing Q particles, all of them 
moving with interactions. To this purpose we make the 
assumption that the density of tracer particles is dilute, 
so that we can consider each tracer particle as though 
it were the only such particle in each cell at any given 
time. This is consistent with processes of self-diffusion 
and mutual diffusion when the tracer is not identical to 
the other particles. The issue of considering many tagged 
particles with correlated motions, such as in a process of 
color diffusion will be discussed in the Conclusions. 

Each cell of our system has a finite size domain which is 
delimited by periodic boundary conditions. This domain 
can for instance be a square in two dimensions or a cube 
in three dimensions. The center of mass is taken to be 
at rest so that no net current takes place. For the sake 
of the argument, we will assume that the periodic array 
of cells has a cylindrical structure, extended along one of 
its spatial dimensions and periodic otherwise. 

The dynamics is thus constructed in a way that it is 
periodic from one cell to the other. A particle entering 
one cell through any of its borders has a corresponding 
image which simultaneously leaves that cell through the 
opposite border, entering the corresponding neighboring 
cell. The cells located at the boundaries of the array have 
however a special status since some of their borders are 
not contiguous to a neighboring cell. Particles can there- 
fore leave the system through these borders, and thus be 
absorbed, while other particles are correspondingly being 
injected through the opposite side of the cylinder, enter- 
ing the corresponding cell through the opposite border. 
While this process occurs at constant rate, it can nev- 
ertheless be used to induce a mass current by simply 
tagging one of the particles at a time with prescribed 
rates of injection into the system. This way a molecu- 
lar dynamics simulation of tracer particle diffusion in a 
spatially extended system can be set up, which involves 
a finite number of degrees of freedom only. 

Thus consider the cubic box of side 1,1 = [—1/2, 1/2Y C 

, containing a collection of Q non-overlapping hard 
sphere particles of radii a with phase-space coordinate 
Tq = {{q%,Pi)}i<i<Q, i e I, e M'*. Let C„ denote 
the phase space of the nth elementary cell. The local 
dynamics, denoted c/j^Tq, preserves the micro-canonical 
measure and is assumed to be chaotic and mixing. 

Given a collection of N copies of such a box, which 
are placed side by side along the first coordinate axis, 
the position of a given tracer particle is specified by its 
coordinate q = {q^ , . . . , q'^) within the elementary box, 
identical to one out of the Q particles there, and by its 
lattice position n, —N/2 < n < N/2 (we assume N even), 
which distinguishes it from its images in the other lattice 
cells. Its real coordinate is therefore {q^+nl, . . . , q"^). The 
flow Fq} acts over the extended system U„C„, i.e. 

on the NQ particles of the N elementary cells, out of 
which a single particle -say the first among Q- has been 
tagged, the corresponding image being located in cell n, 
until it moves to a neighboring cell under the time evo- 
lution $* . Figure [5] shows an example of such a system 
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where the number of tracer particles has however been 
exaggerated for visual purposes. We are indeed assum- 
ing a dilute limit according to which no correlations take 
place between tracer particles. 



FIG. 6: Example of a spatially periodic system with many 
interacting particles and one tagged particle in each cell. The 
tagged particles are randomly injected at the left and right 
boundaries with suitable probabilities p±. For purposes of 
visualization, the number of tracer particles has been exag- 
gerated, mimicking a dilute system where different tracers 
move independently of one another. 



Given tracer injection rates and at the left and right borders, the stationary tracer density is written in a 
way similar to Eq. ()10|1 . in the form 



P{{n,^Q}) = 



P++ P- , P+ - P- 



OO 

i({",rQ}) + ^ [L($-*4n,rQ}) - L{^-'^-^{n,rQ}) 



k=l 



(26) 



This density is assumed to be normalized with respect 
to the local micro-canonical distribution of Q particles 
per cell, and the tracer position along the cylinder axis is 
L({n,rQ}) = nl + q^. We can take the discrete position 
L instead of the continuous position along the horizontal 
axis because both give equivalent results since the frac- 
tal structures appear through the continuous phase-space 
dependence of p on Tq . The factorization occurs because 
the tracer dynamics is a simple passive advection process 
under the assumption that tracer particles do not inter- 
act. Thus equilibration of the interacting particles takes 
place regardless of the tracer dynamics. 

The tracer density ([26]) is a conditional distribution of 
finding a single tagged particle at the phase point Tq 
in cell n given the relative injection rates at the sys- 
tem's boundaries. This density is defined with respect 
to the micro-canonical equilibrium distribution in every 
cell. The continuum limit is recovered, as with the open 
Lorentz gas and multi-baker map, by considering an infi- 
nite number of copies of independent systems, each with 
a single tagged particle. The average number of tagged 
particles as a function of the position along the posi- 
tion coordinate thus yields a macroscopic density which 
evolves according to the Fokker-Planck equation ([3|) . The 
tracer dynamics we consider here is therefore much like 
a single particle system, except for the dimensionality 
of the phase-space dynamics. In particular the diffusion 
coefficient depends on the collective motion of the Q in- 
teracting particles per periodic cell. 

A more physically consistent limit would be to consider 
the absence of correlations between tracer particles set 



in motion in the same arbitrary large system. However 
our formalism for the identification of fractal structures 
in the non-equilibrium stationary states and the compu- 
tation of the entropy production rate makes systematic 
use of the strict spatial periodicity of the dynamics. Its 
extension to such more general situation requires appro- 
priate elaborations. We leave it as a perspective for the 
time being. 

VI. ENTROPY PRODUCTION 

Though a fine-grained Gibbs-type entropy associated 
to a time-dependent phase-space density p(r,t), 

5*(c„) = -/ dTp{r,t)[\ogp{r,t)-i], (27) 

is preserved by the time evolution, this is not the case 
for coarse-grained entropies. This observation has long 
been understood, see e.g. [H, UHl, however the novelty 
here is that, even though the system under consideration 
is Hamiltonian, the stationary state distribution has a 
singular density. The fractal structure of the stationary 
state forbids the use of an entropy like (P?]) . Accordingly, 
the proper treatment of these phase-space measures re- 
quires the use of coarse graining methods. 

A systematic approach to defining the proper coarse- 
grained entropy was outlined in [j, [2^ . The idea is 
that, in the stationary state, owing to the singularity 
of the invariant density, the Gibbs entropy (P7|) should 
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be evaluated by using a grid of phase space, or partition, 
G = {dTj}, into small volume elements dTj, and a time- 
dependent state fin {dTj , t) . The entropy associated to 
cell C„, coarse grained with respect that grid, is defined 
according to 



^fln{drj,t) 



log 



fi„{dTj,t) 
dTi 



1 



(28) 



This entropy changes in a time interval t according to 
A-5*(C„) - 5S(C„)-5^-^(C„), (29) 

where, in the second line, the collection of partition ele- 
ments {dTj} was mapped to {^'^dVj}, which forms a par- 
tition $'^G whose elements are typically stretched along 
the unstable foliations and folded along the stable folia- 
tions. 

Following 0, and in a way analogous to the phe- 
nomenological approach to entropy production fld^ . the 
rate of entropy change can be further decomposed into 
entropy flux and production terms according to 

A-SUCn) = AlSliCn) + AlSliCn) , (30) 

where the entropy flux is defined as the difference be- 
tween the entropy that enters cell C„ and the entropy 
that exits that cell, 

A:4(C„) = 5|$.,r^.j(C„) - 5*^.,r^j($-C„) . (31) 

Collecting Eqs. ([^ - (j^ . the entropy production rate at 
C„ measured with respect to the partition G, is identified 
as 



(32) 



This formula is equally valid in the non-equilibrium sta- 
tionary state. 

This ab initio derivation of the entropy production 
rate for non-equilibrium systems with steady mass cur- 
rents yields results which are in agreement with the phe- 
nomenological prescription of thermodynamics. Indeed, 
as we show below, the entropy production rate of cell 
C„, Eq. ([5^ . reduces to the phenomenological entropy 
production at corresponding position Xn, Eq. ([?]), as the 
grid elements become small and the dependence on the 
choice of partition thus disappears. 



A. Multi-baker map 

This formalism is particularly transparent for the 
multi-baker map. Referring to Fig.Ul and having in mind 
that measures are time-evolved by the inverse map Bq^, 
it is easy to see that the bottom horizontal half of cell n is 
mapped to the left vertical half of cell n+1 and, likewise, 
the top horizontal half of cell n is mapped to the right 



vertical half of cell n — 1. Now the invariant measure 
(|24p is uniform with respect to the x coordinate, which 
is the expanding direction, and has a fractal part along 
the y axis, the contracting direction. Thus the entropy 
of the vertical half of cell n + 1 is half of the entropy of 
cell n + I, but is not equal to the entropy of the bottom 
or top horizontal halves of that cell, at least not with 
respect to the same partition. 

This observation yields the identification of the entropy 
production rate, namely the difference between the en- 
tropies of that cell, measured at two different levels of 
resolution, one with respect to volumes dV — dxdy, the 
other with respect to volumes dV = dx'dy' , where dy' is 
twice as large as dy, and dx' half as large as dx. 

To be precise, given a resolution level 2~*'', k E N, 
we consider the collection of cylinder sets dTk{yj) — 
{{x,y)\yj < y < yj+ 2-^1}, j = l,...,2^ with = 
2 '^(j — 1)1, which partition the square into 2*^ horizon- 
tal slabs of widths 2^'^l, and, by analogy with Eq. 
define the fc-entropy of cell n to be 



5'fe(C„) = - y^jinidTkiyj)) 



log 



(33) 

According to Eq. (|32[) , the fc-entropy production rate per 
unit time is 



AlSkiCn) 



1 

T 

2*3 + 1 



[5fe(C„)-^fe+i(C„)], (34) 

= V t^n{drk+i{yj)) log — — , 

f^n{dTk{yj/2)) 

where the index j/2 is the largest integer above j/2. 

The thermodynamic entropy production rate is recov- 
ered in the continuum limit, whereby the local gradients 
tend to zero and fc can be arbitrarily large. 

To verify this, we consider the symbolic sequences 
iVf. = {^0, . . . , Wfc-i}, i^Jj e {0,1} associated to the 2'^ 
points yj through the dyadic expansions yj = y{^k) — 

In terms of the symbolic sequences w^, the invari- 
ant measure associated to the cylinder set dY{uij.) = 
dT{y{uif.)) becomes, using Eq. 



Mri(Wfc) 



2"V. 



VpAT(^, 



1 



-^2'^AT(^fc) 



(35) 



where we wrote the local density gradient Vp = (p+ — 
p.)l/L, and AT(o.J = T(y(^,)// + 2"'=) - T{y{u^)/l). 
Using the functional equation (|25p . it is straightforward 
to check that 



2'=AT(a;,) = ^(l-2c.,), 

3=0 



(36) 



which, up to the scale I, is the displacement associated 
to points in the cylinder set dT{uj_^,). Two identities are 
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immediate 



2'=AT(a;,) = 0, 



(37) 



Substituting this expression into Eq. ([34]) . we arrive, 
after expanding in powers of the density gradient, to the 
expression 



1 (Vp)2 



2t fin 

1 (Vp)2 ^2 (p^_p_)2 



2t 2tL^ 



(38) 



Identifying the diffusion coefficient T) = PD/t, with 
D = 1/2 the diffusion coefficient associated to the bi- 
nary random walk, we recover the phenomenological en- 
tropy production rate for the macroscopic position 
variable X = nl/L : 



lim Ar5fc(C„) = 



djSjXn = nl) 
dt 



(39) 



where the limit t, Z — > is take with the ration P /t fixed. 



B. Lorentz gas 

This computation transposes verbatim to the non- 
equilibrium stationary state of the Lorentz gas, Eq. (fTO|) . 
whose entropy production rate is computed as follows. 

Let M be a partition of 9C, the phase space of the 
collision map, into non-overlapping sets Aj, M = ^jAj. 
We assume that each element of the partition is contained 
in a cell 9C„, i.e. VjBn : Aj C 9C„. There is thus a 
sub-collection M„ of cells Aj which partition (?C„ , 9C„ = 
UjiA eM„^j- Moreover we assume all the partitions M„ 
are isomorphic, which is to say they can be obtained from 
one another by translation. 

We construct such a partition of the phase space as- 
sociated to cell n by dividing up i9C„ into sets Aj ac- 
cording to the displacements of points F g Aj, such as 
shown in Fig. [T] Thus a k partition M.^ of C„ is a col- 
lection of sets A(uJo, . . . jUJk-i) of points F which have 
coherent horizontal displacements for the first k steps. 
These displacements are coded by sequences of symbols 
LU/. which, in this case, take at most 12 possible values, 
corresponding to the twelve possible transitions to the 
nearest and next-nearest neighboring disks. However, 
unlike with Bernoulli processes, not all symbol sequences 
are allowed so that the number of sets in the partition 
is much less than 12^^. This is referred to as pruning 

[13, m. 

Assuming G {1, . . . , 12}, we define a{ujj) as the dis- 
placement along the horizontal axis corresponding to the 



^ 0- 




FIG. 7: The partition of a unit cell of the open periodic 
Lorentz gas constructed upon two iterations of the collision 
map. The borders of the elements of the partition are the 
lines of discontinuity of the collision map in the {(t>,^) plane. 
Each point of the curves belongs to a trajectory that grazes 
a disk in one or two iterates. 



label ujj , and 



k~l 

^a(wj). 



(40) 



which, up to the local length scale /, is the displacement 
associated to the sequence w , measured with respect to 
the position of the scatterers. 

The measure associated to these sets is defined by 



fiiyj,) = / dFp(F) 



(41) 



where p is the invariant density associated to the fiux 
boundary conditions, Eq. (jlOp . Inserting this expression 
into Eq. (j4ip . we may write 



M(^fc) = ^J■nv{oJ.k) (42) 

p CO 

_^P±_P_ / dVy][x{^-'''T) - x{^-'^-^V)] 
L Ja(ui,:) 

where 

P++ P- 



fe=i 



Pn = 



iP+ - P-)^ 



(43) 



and v{oJ_i?) — J^^^ ^ dT is the volume measure associated 
to the set A{uj_,^), with i'(wj,) = whenever w^, is a pruned 
sequence. 

We claim that we can make the approximation 



Siulk) , (44) 
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which we expect to hold provided k is large. This quan- 
tity therefore plays, for the non-equilibrium stationary 
state of the Lorentz g role similar to that played un- 
der the same conditions by AT{uj_f.)/i^{uj_f.), Eq. ([55]) . for 
the multi-baker map. Indeed the difference between the 
two sides of Eq. is 0(1). Under this condition and 
assuming / ^ 1, we thus write 



1 



(45) 



which can be compared to Eq. (|35p. When k is large, 
identities similar to Eq. (I37p hold for (5(wj.) : 



E.,i^(^J'5(c.,)-o, 

E., K^J<5(^fe)2-2i^fc, 



(46) 



where the diffusion coefficient D is that of the random 
walk on half integer lattice sites associated to the dis- 
placements PO)) with weights i^(wj.). 

The computation of the entropy production rate (j32[) 
then proceeds along the lines of Eq. ([55]) , with a similar 
result: 



Ar5fe(C„ 



(47) 



in agreement with the phenomenological expression 
Eq. ([T3|) with diffusion coefhcient V = PD/t and po- 
sition Xn — nl/2. 



C. Interacting particle system 

Here, as in Ref. 0], we assume a partition of phase- 
space in sets Aj sufficiently small that all points in them 
have particle No. 1 flowing through the same sequence of 
cells over a large time interval, i.e. such that the set of 
points ^~'^Aj are in the same cell when projected along 
the coordinate of particle No. 1, with < t < T. The 
corresponding location is L($^'^{n, Fq}), Tq e Aj. We 
define the displacement 



d{TQ,t)^L{^-*{n,TQ})-nl, 
and write the non-equilibrium stationary state as 



(48) 



p{Aj) ^flnV{Aj) 



p+ - P- 



dVQA{TQ,r), (49) 



where the volume measure v{Aj) is the micro-canonical 
measure of the periodic system of Q particles. The inte- 
gral in this expression encodes the fluctuating part of the 
non-equilibrium stationary state. In analogy to Eq. (|44ll . 
we can write 



^(A,)^^/^dr,d(r„r), 



in terms of which we have the identities 

Y.,y{Aj)5{Aj) = Q., 
Y.,y{A,)5{Ajf^2VT, 



(50) 



(51) 



which, apart for their dimensions, are similar to Eqs. (|37[) 
and 

Proceeding as with the multi-baker map and Lorentz 
gas, we have the entropy production rate 



Ar^T(C„) 



P (P+ -P-) 



(52) 



independently of the resolution scale set by the time pa- 
rameter T. 



VII. CONCLUSION 

In this paper we have considered a large class of 
diffusive deterministic systems with volume-preserving 
dynamics, and established the fractality of the non- 
equilibrium stationary states which result from the 
stochastic injection of particles at the systems' bound- 
aries. 

Under the assumption that the dynamics is chaotic, 
the natural non-equilibrium invariant measure is charac- 
terized by a phase-space density which is the stationary 
state of a Liouville or Perron-Frobenius operator with 
flux boundary conditions, which account for the pres- 
ence of particle reservoirs. This density naturally splits 
into regular and singular parts. The regular part has 
a local equilibrium form, uniform with respect to the 
scale of macroscopic variables. It is the natural coun- 
terpart on phase space of the solution of the macroscopic 
transport equation under non-equilibrium boundary con- 
ditions, displaying a linear density gradient. In contrast, 
the singular part varies non-continuously on microscopic 
scales, reflecting the sensitive dependence on initial con- 
ditions, and has no macroscopic counterpart. 

The fundamental achievement of this paper is to have 
exhibited the universal features of the non-equilibrium 
stationary states of simple low-dimensional deterministic 
models of diffusion, which are shared by higher dimen- 
sional spatially periodic systems modeling tagged particle 
diffusion. As already pointed out in ^T3| , the fractality of 
the non-equilibrium stationary states is the key to a sys- 
tematic computation of the positive entropy production 
rate. Indeed the fractality of phase-space distributions 
imposes the use of coarse-graining techniques in order to 
define the entropy. The coarse-graining of phase-space 
into partitions of sets of strictly positive volumes results 
into an entropy which is generally not constant in time. 
The decomposition of the time evolution of such coarse- 
grained entropies into entropy flux and entropy produc- 
tion terms provides a framework in which the positiveness 
of the entropy production rate can be rigorously estab- 
lished. Furthermore, in the macroscopic limit, two key 
results are obtained: (i) the entropy production rate is 
independent of the coarse graining in the limit of arbi- 
trarily flne phase-space graining, and (ii) its value is iden- 
tical to the phenomenological rate of entropy production, 
according to thermodynamics. 



13 



In this respect, this paper helps clarify a point that 
had been overlooked in previous works. Under the diffu- 
sive scaling limit, the local particle density gradient is a 
natural small parameter in the continuous limit: fixing 
the injection rates and system size, wc let the spacing be- 
tween cells go to zero in order to obtain the macroscopic 
limit. With the identification of this small parameter, 
it is straightforward to verify that the computation of 
the entropy production of the non-equilibrium stationary 
states, whether of multi-baker maps, Lorentz gases, or 
spatially periodic systems, all yield leading contributions 
consistent with the prescription of thermodynamics. 

As we have pointed out, the formalism we use relies 
systematically on the spatial periodicity of the dynamics, 
especially with regards to identifying the rate of entropy 
production. Though this periodicity is convenient as it 
enables one to clearly identify a separation of scales be- 
tween microscopic and macroscopic motions, it is a rather 
restrictive and ultimately undesirable feature of our mod- 
els. One would rather describe the diffusive motion of di- 
lute tracer particles in arbitrary systems. It is our hope 
that this work helps set the stage to achieve this more 
ambitious goal in future works. 

As mentioned in the introduction, this paper is the first 
of two. In the companion paper [12| . we consider the in- 



fluence of an external field on spatially periodic diffusive 
volume-preserving systems, i.e. forced systems in which 
no dissipative mechanism is present. As proved by Cher- 
nov and Dolgopyat [l5|, such systems are recurrent in 
the sense that tracer particles keep coming back to the 
region of near zero velocity where all the energy is trans- 
ferred into potential energy, so that no net current takes 
place. We show that these systems remain diffusive, al- 
beit with a velocity-dependent diffusion coefficient which 
alters the usual scaling laws and consider their statistical 
properties in some details. 
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